///@file ProbUtils.cpp
///@brief Impementation of utility functions for statistical studies
///@author Arnaud Duval
///@version 0.0
///@date 2010/10/27

#include <assert.h>
#include <cmath>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

#include "ProbUtils.h"

///@brief Compute factorial of n
///@param n Number for which the factorial must be computed
///@warning This function is limited to n=12
///@return Value of n!
unsigned long int Factorial(int n)
{
	assert(n < 13);
	unsigned long int res = 1;
	for(unsigned long int i = n ; i > 0 ; i--)
		res *= i;
	return res;
}


///@brief Compute binomial coefficient
///@param n Cardinality of the set
///@param k Cardinality of the selected set
///@return Binomial coefficient C_n^k
unsigned long int Binomial(int n, int k)
{
	assert(k <= n);
	if(k == 0)
		return 1u;
	else if(k == n)
		return 1u;
	else
		return Binomial(n-1, k)+Binomial(n-1, k-1);
}


///@brief Compute an approximation of reciprocal of erf functions
///@param z value for which the reciprocal should be computed
///@param o order of approximation (a value of 1000 gives good results)
///@return Value of erf reciprocal
double Recerf(double z, int o)
{
	double res = 0.;
	double *c = new double[o+1];
	c[0] = 1.;
	for(int k = 0 ; k <= o ; k++)
	{
		if(k != 0)
		{
			c[k] = 0.;
			for(int m = 0 ; m<= k-1 ; m++)
			{
				c[k] += (c[m]*c[k-1-m])/((m+1.)*(2.*m + 1.));
			}
		}
		res += (c[k]/(2.*k + 1.))*pow((sqrt(M_PI)*z)/2.,2.*k + 1.);
	}

	delete[] c;
	return res;
}

///@brief Compute distribution function of a uniform random variable
///@param x Value for which the function should be computed
///@param a lower bound
///@param b upper bound
///@return Value of the distribution of a uniform random variable
double Fu(double x, double a, double b)
{
	assert(a < b);
	if((x >= a)&&(x <= b))
		return (x-a)/(b-a);
	else
		return 0.;
}

///@brief Compute reciprocal of uniform distribution
///@param x Value for which the function should be computed
///@param a lower bound
///@param b upper bound
///@return Value of the reciprocal of a uniform distribution
double Recu(double x, double a, double b)
{
	assert(a < b);
	return (b-a)*x + a;
}

///@brief Compute reciprocal of Gauss distribution function
///@param x Value for which the function should be computed
///@param mu average value of Gauss distribution
///@param sigma standard deviation of Gauss distribution
double RecFgauss(double x, double mu, double sigma)
{
	return (sigma * sqrt(2.) * Recerf(2.*x - 1., 1000)) + mu;
}
